Quantum Monte Carlo simulation of overpressurized liquid ^He 
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A difFusion Monte Carlo simulation of superfluid ''He at zero temperature and pressures up to 
275 bar is presented. Increasing the pressure beyond freezing (~ 25 bar), the liquid enters the 
overpressurized phase in a metastable state. In this regime, we report results of the equation of 
state and the pressure dependence of the static structure factor, the condensate fraction, and the 
excited-state energy corresponding to the roton. Along this large pressure range, both the condensate 
fraction and the roton energy decrease but do not become zero. The roton energies obtained are 
compared with recent experimental data in the overpressurized regime. 

PACS numbers: 67.40.-w,67.80.-s 



Quantum fluids in metastable states are presently a 
research topic of fundamental interest from both the ex- 
perimental and theoretical viewpoints . The extremely 
low temperature achieved in liquid helium makes this liq- 
uid the most pure in nature and therefore the optimal 
choice for observing homogeneous nucleation, which is 
an intrinsic property of the liquid. Caupin, Balibar and 
collaborators have studied profusely the negative pres- 
sure regime by focusing high intensity ultrasound bursts 
in bulk helium liquid "^He they have measured a 

negative pressure of -9.4 bar, only 0.2 bar above the spin- 
odal point predicted by microscopic theory 0, Q . The 
same experimental team has used recently this acoustic 
technique to pressurize small quantities of liquid ^He up 
to 160 bar at temperatures 0.05 K < T < 1 K 0. This 
pressure is the highest pressure ever realized in overpres- 
surized liquid '*He and is much larger than the liquid- 
solid equilibrium pressure, which atT = 0Kis25.3 bar. 
In the experimental setup used by Werner et al. 6j, the 
overpressurized regime has become accessible by avoid- 
ing nucleation on the walls of the container and therefore 
allowing only homogeneous nucleation. This nucleation 
was not observed along this large increase of the liquid 
pressure beyond the freezing point. 

Liquid ^He in metastable states has also been obtained 
by immersing it in different porous media. Albergamo et 
al. have carried out neutron scattering experiments 
in a medium with 47 A pore diameter filled with den- 
sities below the equilibrium density and negative pres- 
sures up to -5 bar. The achievement of negative pres- 
sures in this medium is attributed to the stretching that 
the liquid supports due to the strong attraction to the 
walls. Using a different material, with 44 A pore diam- 
eter, Pearce et al. have reported neutron scattering 
data in the high density regime observing a liquid phase 
up to ^40 bar. Therefore, the confinement of helium in 
porous media makes feasible extensions of the pressure 



on both sides of the stable liquid phase. The nature of 
these two metastable regions presents a significant dif- 
ference. At negative pressure, there exists an end point 
(spinodal point) where the speed of sound becomes zero, 
and it is thermodynamically forbidden to cross it main- 
taining a homogeneous liquid phase. On the contrary, 
such a point does not exist on the overpressurized side. 
Nevertheless, Schneider and Enz suggested that the 
pressurized phase has also an end point corresponding 
to the pressure where the excitation energy of the roton 
might vanish. On the other hand, and according to Jack- 
son et al. ^■iid Halinen et al. i the vanishing of the 
roton energy is not the instability that causes the liquid- 
solid transition in "^He but another one involving a 6-fold 
symmetric soft mode in the two-body correlations. 

The roton mode is not observed anymore when liquid 
"^He crystallizes. In the overpressurized regime, between 
25 and 38.5 bar, the neutron-scattering data of Ref. Q 
show the roton excitation, while the maxon disappears, 
probably due to the fact that its corresponding energy 
might exceed twice the roton energy. At the liquid-solid 
transition, the roton energy is still finite. The latter fea- 
ture means that *He is superfluid when it crystallizes. 
This result seems to be in disagreement with a recent ex- 
perimental work by Yamamoto et al. [l^ who reported 
a critical temperature Tc = at a pressure P ~ 35 bar 
in a porous material. As a first explanation, the differ- 
ences observed in the two experiments can be attributed 
to the smaller pore diameter (25 A) used by Yamamoto 
et al. J but additional work would be necessary to con- 
firm this argument. It is worth noticing that the findings 
of Ref. J2] imply the existence of a superfiuid to nor- 
mal phase transition at zero temperature, an intriguing 
possibility that makes the overpressurized regime even 
more interesting. In this direction, Nozieres 13] has pre- 
dicted recently that the condensate fraction could vanish 
at a certain pressure, and therefore a normal liquid before 
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freezing would be possible. 

At present, the theoretical knowledge of the metastable 
regime at negative pressures is rather complete with an 
overall agreement among quantum Monte Carlo, hyper- 
netted chain (HNC) based on Euler-Lagrange optimiza- 
tion, and density functional approaches [l| . On the con- 
trary, the pressurized liquid remains until now nearly un- 
explored. In the present work, we present diffusion Monte 
Carlo (DMC) results for the overpressurized phase up to 
P ~ 275 bar. The DMC method is probably the best 
suited way to deal with this metastable regime since the 
physical phase of the system is controlled by the trial 
wave function used for importance sampling. Our re- 
sults show a superfluid phase in which, up to the highest 
pressure studied, the condensate fraction is small but dis- 
cernible and the excitation energy of the roton decreases 
but does not reach the zero limit. 

The DMC method is nowadays a well-known tool de- 
vised to study quantum fluids and solids at zero tem- 
perature. Its starting point is the Schrodinger equation 
written in imaginary time. 



' dt 



= {H-E,)^{R,t) 



with an A'^-particle Hamiltonian 
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In Eq. JQ, E^. is a constant acting as a reference en- 
ergy and R = (ri, . . . ,rff) is a walker in Monte Carlo 
terminology. DMC solves stochastically the Schrodinger 
equation ^ replacing '^{R, t) by t) = ^{R, t)ip{R) 
with tp{R) a trial wave function used for importance sam- 
pling. When i — » oo only the lowest energy eigenfunc- 
tion, not orthogonal to 4^{R), survives. The simulation 
of the liquid in its ground state is carried out by using 
a Jastrow approximation, ipjiR) — Y[i<j fi''~ij)- 
our previous work 0| on the equation of state of liq- 
uid ''He in the stable regime, /(r) is a model proposed 
by Reatto [lH which incorporates nearly optimal short 
and medium range two-body correlations. As a matter of 
comparison, simulations of the crystalline hep phase have 
been also carried out; in this case, a Nosanow- Jastrow 
model is used, ipN.i{R) = ij.j{R)Yl^ g{ni), with g{r) a 
gaussian function linking every particle i to a fixed point 
n of the lattice. 

The calculation of the roton energy using the DMC 
method is more involved since it corresponds to an ex- 
cited state. In this case, we have taken as a trial wave 
function for importance sampling an eigenstate of the 
total momentum operator which incorporates backflow 
correlations, as originally proposed by Feynman and Co- 
hen 
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FIG. 1: Energy per particle of liquid ^He from the equilibrium 
density up to the highest density calculated, 0.6 a~^, (solid 
circles). The solid line corresponds to the fit to the DMC 
energies using Eq. and the open circles are experimental 
data in the stable regime from Ref. DMC results for 

the solid phase are shown as solid squares and compared with 
experimental data from Ref. fTsIl (open squares). The error 
bars of our data are smaller than the size of the symbols. 



with fi = ri + Y.j'^iV{ri])rij, and 7y(r) = Aexp[-((r - 
rb)/wb)^]- The only objective of our calculation is the 
energy of the collective mode, which is the same for ex- 
citations with momenta q and —q. Therefore, we avoid 
the complexity of working with a complex wave function 
© by using a superposition of both states. Proceeding 
in this form, the excited state turns into a fermion-like 
problem since the trial wave function is real but not pos- 
itive everywhere ^ fi^'^t step, we have used the 
fixed-node (FN) approximation, which provides an upper 
bound to the roton energy. We verified that the intro- 
duction of backfiow correlations in the trial wave function 
produces results quite close to experimental data at the 
equilibrium density, especially near the roton minimum. 
The nodal constraint imposed by FN is removed, in a sec- 
ond step, by using the released-node (RN) technique in 
which the walkers are allowed to cross the nodal surface 
imposed by tp a-nd survive for a finite lifetime. 

In Fig. 1, we report the results obtained for the en- 
ergy per particle of the liquid phase as a function of the 
density. The He-He interaction corresponds to the HFD- 
B(HE) Aziz potential and the number of atoms used 
in the simulation increases with density: A'^ = 150 near 
freezing and A^ = 250 at the highest densities; periodic 
boundary conditions are assumed in all the simulations. 
Possible biases coming from the use of a finite number 
of walkers and from a finite time-step in the employed 
second-order algorithm are reduced to the level of the 
statistical noise. The DMC energies are accurately pa- 
rameterized, from the spinodal point up to the highest 
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FIG. 2: Pressure as a function of the density. The sohd 
hues stand for the DMC results obtained from the equa- 
tions of state of the liquid and solid phases shown in Fig. 
1. The dashed line is the extrapolation from experimental 
data 6, 22] ; the symbols correspond to experimental data for 
the liquid 1^ and solid phases 

densities calculated, by the analytical form 

e{p) = eo + ei{p/pc-l){l- {p/pc-l)) 

+b,{p/p,-lf + b4p/p,-l)\ (4) 

with e = E/N, and pc = 0.264 a'^ (cr = 2.556 A) the 
spinodal density. The rest of parameters in Eq. Q are 
eo = -6.3884(40) K, ei = -4.274(31) K, 63 = 1.532(12) 
K, and 64 = 1.433(24) the figures in parenthesis being 
the statistical errors \M DMC results of the energies 
of the solid phase, calculated using the Nosanow- Jastrow 
trial wave function and an hep lattice, are also shown 
in Fig. 1. The comparison between the liquid and solid 
phase simulations shows clearly that DMC is effectively 
able to study the overpressurized liquid phase in spite of 
not being the ground-state (minimum energy) configu- 
ration, which manifestly corresponds to the solid phase 
beyond the freezing point. 

Using the equation of state Q), the pressure is obtained 
from its thermodynamic definition P{p) — p^{de/ dp). 
The results obtained, which are shown in Fig. 2, repro- 
duce accurately the experimental data in the stable 
regime and predict a pressure P ~ 275 bar at the highest 
density evaluated, p = 0.6 . Our results are compared 
in the same figure with the analytic form suggested in 
Ref. 0, adjusted to Abraham's experimental data [2^ . 
Below the freezing point, both curves agree but they give 
significantly different values at higher densities; the dif- 
ference amounts to ^ 100 bar at p = 0.6 . As a 
matter of comparison, the figure also shows the pressure 
of the solid phase, derived from the DMC equation of 
state shown in Fig. 1. 

A characteristic feature of a solid phase is the presence 
of high-intensity peaks in the static structure function 
S{q) = {pqp^q)/N, pq = X^ili t^c rcciprocal 

lattice sites. Following the overpressurized liquid phase. 
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FIG. 3: Static structure function of the liquid phase for differ- 
ent densities. From bottom to top in the height of the main 
peak, the results correspond to densities 0.365, 0.438, 0.490, 
0.540, and 0.6 a'^. 
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FIG. 4: Condensate fraction of liquid ''He in the overpres- 
surized regime. The line is an exponential fit to the DMC 
results. 

we have not observed this feature and thus the liquid na- 
ture of the system is confirmed. In Fig. 3, results of S{q) 
for densities ranging from equilibrium up to the highest 
densities studied are reported. The results show the ex- 
pected behavior: when p increases, the strength of the 
main peak increases and moves to higher momenta in a 
monotonic way. At low momenta, the slope of S{q) de- 
creases with the density, following the limiting behavior 
limq^o S{q) = hq/{2mc) driven by the speed of sound c. 

A characteristic signature of bulk superfluid ^He is a 
finite value of its condensate fraction, i.e., the fraction of 
particles occupying the zero-momentum state. As usual 
in a homogeneous system, the condensate fraction no has 
been estimated from the long range behavior of the one- 
body density matrix, limr_»oo p('') — "-o- The results 
obtained for no, from the melting pressure up to nearly 
300 bar, are plotted in Fig. 4. The line on top of the 
data corresponds to an exponential fit which reproduces 
quite accurately our DMC results (no(P) ~ Ae'^^, with 
A = 0.052 and /3 = 0.015 bar"-'^). As one can see in 
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FIG. 5: Roton energy as a function of the density (solid cir- 
cles). Open circles stand for experimental data from Ref. 0. 
The line is an exponential fit to the DMC data. 



the figure, no decreases quite fast until P = 100 bar and 
then the slope decreases, approaching a value uq ~ 0.005 
at the highest density. With the same procedure, we 
obtained 14] uq — 0.084(1) at the equilibrium density 
Pq — 0.365 <J~^, value which is compatible with PIMC 
estimations at low temperature 0| (0.069(10) at T = 
1.18 K and 0.087(10) at T = 1.54 K). 

Superfluidity in bulk "^He is manifested in the dynamic 
structure function by the clear signature of the roton col- 
lective mode. By increasing the density the roton energy 
Cr decreases, and its approach to zero has been proposed 
as a possible final point in the overpressurized regime . 
Trying to discern this hypothesis, we have carried out a 
RN calculation of the roton energy beyond the freezing 
point. The same methodology was used in the past in a 
DMC calculation of the phonon-roton spectrum at equi- 
librium and freezing densities [T^ arriving at an accurate 
description of the experimental data. The results ob- 
tained are shown in Fig. 5. Experimental data obtained 
by neutron scattering experiments on superfluid ^He in 
a porous media and up to 40 bar are also plotted in the 
figure p. From to 40 bar, decreases linearly with 
the pressure and our data reproduces well this behavior. 
However, increasing the pressure this slope is reduced 
and, at the highest density studied the roton energy is 
still different from zero (e,. = 2.8±1.2 K at p = 0.58 a~^). 
It is worth noticing that the statistical errors are rather 
large and difficult to reduce since er = E,- — Eq, with 
JEJr and Eg the total energy of the excited and ground 
state, respectively. The excited state energy is estimated 
through an exponential fit E{t) — E^ + Ae~^*'^^\ with t 
the released time. The uncertainty of this extrapolation 
is under control since in all cases the difference between 
considering the last calculated point in released time or 
Ef is of the same order as the statistical noise. At each 
density, the number of particles has been adjusted to be 
as close to the roton momentum as possible; neverthe- 



less, as only discrete values of q are accessible, possible 
corrections to the energy due to this fact are estimated to 
be less than 0.5 K in all cases. The obtained values of the 
roton momentum do not follow linear dependence with 
pressure, as the data in the stable liquid regime might 
suggest |2J. As the pressure is increased, the slope of 
the roton momentum as a function of pressure decreases. 

In summary, the present DMC study of superfluid ''He 
in the overpressurized regime does not show any sig- 
nature of final (spinodal) point. The structure factor, 
condensate fraction, and roton energy are mainly driven 
by the density and not by the pressure. As shown in 
Figs. 1 and 2, equal changes in density in the stable and 
metastable regimes produce different changes in pressure. 
This leads, in the density range studied, to an approxi- 
mated exponential behavior with pressure of magnitudes 
like no and e^- 
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